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ABSTRACT 



One case of analytic data and one case of real data were numerically 
integrated using a 5~level baroclinic primitive equation global model of 
the general circulation* The feasibility of using a mixed second order 
and fourth order space differencing scheme to improve the phase speeds 
of meteorological waves was examined. The results indicate that mixed 
scheme tends to give a better representation of the phase speeds than the 
second order scheme does. 
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I, INTRODUCTION 



A baroclinic primitive equation model using a global, staggered, 
spherical, sigma coordinate system which was originally written by 
Dr, F, J, Winninghoff and which is now being developed by the United 
States Navy Fleet Numerical Weather Central (Elias, 1973) was used in 
this study. The purpose of this investigation was to modify the hori- 
zontal space difference equations to a mixed second and fourth order 
scheme with the aim of improving the phase speeds of the meteorologi- 
cal waves, as proposed by Williams (1972). Real time data (Elias, 
1973) and analytic data using an analytic spherical harmonic stream 
function (Neatam, 1946) were used as initial conditions. The use of 
the analytic case allowed the control of temperature and moisture dis- 
tributions, predominant wave number, phase speed and wave amplitude. 
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II. BARQCLINIC PRIMITIVE EQUATION MODEL 



The governing differential equations, written in vector form, are 
similar to sets used by Smagorinsky et al (1965), Arakawa, et al (1969) 
and Kesel and Winninghoff (1972). The integrations were carried out on 
a global, spherical, staggered, sigma coordinate system using the con- 
servative-type difference equations based on the procedure of Arakawa 
(1966). The complete set of finite difference equations are given in 
Kesel and Winninghoff (1972) for the non-staggered, polar stereographic 
grid. Appendix A contains the finite difference equations for the 
spherical, sigma, staggered system of equations. The heating, moisture, 
friction and diffusion finite difference terms are given in Kesel and 
Winninghoff (1972). 

A. PRIMITIVE EQUATIONS 

The following are the equations for the model in spherical, sigma 
(a) coordinates. 

The dimensionless vertical coordinate is: 



a = P/tt . 



The vertical velocity measure is defined as: 



w 



da 

dt 



a . 



The latitudinal component of the equation of motion is: 



9u7T 

3t 



1 ,3(uu7T) 3(uv7TCOS0) 

a cos 0 ^ 3X 30 
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, 7 t 3 (wu) . TTuvtanS . - 

+ 1^- ^ — T— * 



« l"lf 1x1 ■*■ " ^2 “ 

a cos 0 3 a 3a 2 u 



( 1 ) 



The longitudinal-component of the equation of motion is: 



3 v 7T _ 1 p 3(uV7T) 3 (w 7 TCOs 9) , 

3t a cos 9 3X 39 



, 7 t 3 (wv) - 7Tuutan0 - 

+ H Tfuf 

3o a 



1 .TT 3$ . RT 37T, . ^ ^ 

- 1 ‘ se * w' “v • 



( 2 ) 



The First Law is: 



9ttT _ 1 , 9(ttuT) 9(TrvTcos9) 

9t a cos 0 ^ 9X 96 



Tr9(wT) 

9o~ 



, RT r . -9 tt 

+ y t-“" + “(Tr 



1 fU 97T , 97 T cos 9 tm 

q [ -5T + 

a cos 9 9 a 99 



+ ttH + . 



(3) 



The Moisture Equation is: 



9qTT _ 1 , 9(TTuq) 9(TTvqcos6) -■ 

9t a cos 0 ^ 9X 90 



+ "iisai +Q,i + D . 

9o q 



(4) 
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The continuity equation is: 



3tt _ 1 -.3(u7T) 3(v7TCOS0), 7T 3w 

3t " a cos 0 ^ 3A 30 ^ 3a * 



The Hydrostatic Equation is: 



3a "a 



( 6 ) 



In the above equations, u is the zonal wind component; v, the 
meridional wind component; 7T, the terrain pressure, T, the temperature; 
f, the Coriolis parameter; $, the geopotential; q, the specific humidity; 
H, the diabatic heating; F, the frictional stress; D, the lateral diffu- 
sion of momentum; D^, the lateral diffusion of. heat; D^, the lateral 
diffusion of moisture source/sink term; P, the pressure; a, the radius 
of the earth; X, longitude and 0, the latitude. 

Kesel and Winninghoff (1972) modified the pressure gradient terms 
as suggested by Kurihara (1968) to reduce the inconsistent truncation 
error which arises as a result of the sensitivity of the model to varia- 
tions in the height and character of the terrain. The value of the geo- 
potential is interpolated at neighboring gridpoint locations to the same 
pressure surface as at the computation point to enable one to compute 
to geopotential gradients on pressure surfaces which are thus synthe- 
sized locally around each computational point. The pressure gradient 
term in sigma coordinates is: 

_ V $ = - V $ — Vti . (7) 

p an 
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Solving for the gradient of tt in component form one gets, 



^ = JL- - — 1 

9A RT ^9X 9X^'' 

P o 



( 8 ) 



in = IL. r— — 1 

90 RT ^96 ” 99 

P o 



(9) 



It is shown in Appendix A that the right hand side of equations (8) and 
(9) are the sum of the thicknesses between the sigma surface and the 
pressure surface to the right and left of the computational points. 

Since the thickness depends only upon the mean temperature of the layer 
in question, the thickness equation is used to solve equations (8) and 
(9). The mean temperature is obtained by interpolation as indicated in 
Appendix A. 

B. GRID 

The longitudinal and latitudinal grid increments were ten degrees 
for the mass and velocity fields. The geopotential ($), temperature 
(T) , specific humidity (q) , and the vertical-velocity (w = -a) were 
carried at the poles. When $, T, q, w and 7T were needed at a velocity 
point they were taken as the average of the four surrounding mass points. 
Figure 1 represents the staggered grid and location of the variables. 

C. VERTICAL LAYERING 

Figure 2 is the vertical layering as described by Elias (1973). 

The basic variables, as indicated, are carried at the center of each 
layer except that the specific humidity (q) is carried at the lower three 
layers only. The vertical velocity (w) is carried at the interface 
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Figure 1. Location of Variables and Grid. 
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Figure 2. Vertical Layering 
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between layers. The vertical coordinate, sigma (a) was defined by 
Phillips (1957). 

D. TIME DIFFERENCING 

The two, finite difference techniques that were used were the cen- 
tered difference scheme and the Matsuno (Euler-backward) difference 
scheme. The finite difference equation for the centered technique is 



F 



t+1 



= F 



t-1 



+ 2At 




and that for the Matsuno scheme is 



F*^ = F*^ + At 



3t 



+ At If 



Since the centered scheme is not feasible for the first time step and 
also produces a computational mode (Haltiner, 1971), the Matsuno scheme 
was used for the initial time step of each 3-hour integration to reduce 
solution separation. The Matsuno scheme also selectively dampens high 
frequency waves (Haltiner, 1971). The Arakawa technique (Langlois and 
Kwok, 1969) of averaging quantities involved in the longitudinal deri- 
vitives is used. This allowed a time step of ten minutes to be used. 
Without this technique, the Von Neuman linear computational stability 
criterion (Haltiner, 1971) would require a two and one half minute time 
step because of the longitudinal grid distance at 85^ latitude which is 
only 47 kilometers. 
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The lateral diffusion, friction, convective condensation and large 
scale condensation are computed at every time step while heating is com- 
puted every six time steps (Kesel and Winninghoff, 1972). 

E. FOURTH ORDER MIXED SPACE DIFFERENCING 

Williams (1972) found that the stability criterion for fourth order 
unstaggered space differencing required a time step which was 27% smaller 
than that for second order space differencing. He then noted that the 
phase speed of the meteorological wave is independent of the difference 
approximation which is used for the pressure force term. A mixture of 
fourth order differencing in the advection terms and the divergence term 
in the continuity equation and second order differencing in the pressure 
force term was proposed. The stability criterion for this mixed scheme 
was found to require a time step that was only 14% smaller than that of 
the second order scheme. Thus the mixed differencing allows a larger 
time step while giving fourth order accuracy in the phase speed of the 
meteorological wave. Appendix B presents the fourth order differencing 
of the advection terms and the divergence term of the continuity equa- 
tion. A time step of 8 minutes was used for this scheme. 
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III. INITIAL CONDITIONS 



Two cases, one using real data and the other analytic data, were 
used for initialization of both the second order and mixed differencing 
schemes. The real data was produced external to the main program on a 
Northern Hemispheric FNWC 63 x 63 grid. The analytic data with known 
phase propagation was produced following the work of Neamton (1946), 
Gates (1962), Heburn (1972) and Elias (1973). The real data were fields 
analyzed by FIWC objective schemes. 

The main program was used to interpolate the data so that values as 
5® latitude-longitude intersections were available on the Northern Hemi- 
sphere (Elias, 1973). The data was then reflected into the Southern 
Hemisphere. In both cases, the input data consisted of temperature 
analyses for the Northern Hemisphere at 12 constant pressure levels from 
1000 to 50 mb, height analyses at ten of these levels, moisture analyses 
at four levels from the surface to 500 mb, sea level pressure and sea 
surface temperature. Monthly mean surface temperature fields were used 
to derive an albedo field (Dickson and Posey, 1967). In both cases the 
terrain height was set to zero. For the barotropic analytic case heat- 
ing, convective adjustment and friction were turned off. 

If one assumes nondivergent horizontal flow, harmonic wave solutions 
of the complete vorticity equation can be obtained for the sphere 
(Neamton, 1946). These solutions give the velocity of propagation of 
the wave. The solution of the field (Elias, 1973) was found to be 

Ip = A sin (mA-vt)P™ (sin 0) - B a sin 6 + ^ (sin 0) (1) 
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where A, B and C are constants to be determined, a is the radius of the 
earth, m is the hemispheric wave number, v/m is the angular phase speed 
of the wave (radians /second) , denotes the Legendre polynomial and 

represents the associated Legendre function. 

If C = 0 and n = m + 1 (Elias, 1973), equation (1) reduces to 

ip = A sin (mA - vt) (2N!/2^N!)sin 0 (cos 0)^ - B a^ sin 0 (2) 



The constant A is arbitrary and proportional to the wave amplitude. It 
was shown by Haurwitz (1946) that the solution obtained for \p implies 
the existence of a velocity distribution over a sphere such that the 
angular velocity of the westerly current is made up of the sum of three 
terms. The first, B, is constant over the sphere; the second term 
varies along a meridian but is constant around a latitude circle; and 
the third term represents the angular phase speed which is given by 



v/m = B 



n(n+l)~2 

n(n+l) 



2fi 

n(n+l) 



(3) 



where Q is the earth’s rotation. For the analytic experiment A was chosen 

2 “1 -7 “1 

as 1000 m sec , m as 6, n as 7 and B as 8 x 10 sec . The solution of 

equation ( 3) gives an angular velocity of ~9® longitude/day (toward the 

west). Since wave number six was used, integrations were performed for 

only one-sixth of a latitude circle with the remaining five-sixths being 

set equal to the first part. This resulted in a 40% saving of computer 

time for the analytic case only. 

The gcopotential fields were derived by solution of the linear 
balance equation and the temperatures, constant at each pressure surface. 
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I. 



i m 




were consistent with the National Advisory Connnittee for Aeronautics 
(NACA) . Both the geopotential and temperature fields were derived as 
shown by Elias (1973). 

The wind field was produced following the work of Elias (1973) by 
using the finite difference expressions 



u 



1 ^ 
a ^6 



(A) 



V = i M 

a cos 6 AX 



(5) 



where a is the earth’s radius, A0 is the latitudinal distance increment 
and AX is the longitudinal distance increment. Poleward of 25® latitude 
the initial geopotential field, was retained and the \p-field was 

obtained by the solution of the linear balance equation, 

2 2 

V + Vi/; • Vf/f = V (|)^/f (6) 



Equatorward of 25® latitude the i|;-field was determined using, 

ip = <I>^/f 

where f is a mean Coriolis parameter. A new geopotential field ($) 
was obtained using 

2 2 

V ^ = (fV \p + Vip • Vf) . 

At 25® latitude, and (f), were a combination of half the value of the 
poleward case and half of the equatorward case. 
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IV. WAVE ANALYSIS METHOD 



A Fourier series was determined at each five degrees of latitude 
around the latitude circle. With this technique, the phase angles and 
amplitudes of each wave number around a latitude circle can be calcu- 
lated. A Fourier series can be expressed as follows (Heburn, 1972): 



F(x) = A + (A cos mx + B sin mx) 

o m m 



m 



= + y ^ (C^ cos (mx - 6m) ) 



m 



B 



where 



C = — 



m 



m 



'm sin(6m) cos (6m) 



and 



6 = tan ^ (B /A ) 

m mm 



The analytic experiment involved an input stream function of wave 
number six. The values of primary interest were 6^ and which are 
the phase angle and amplitude of wave number six. Other angles and 
amplitudes were extracted, especially in the real data case. 
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V. RESULTS 



Two experiments will be presented in this section. The first 
experiment was performed with analytic fields and the second with FNWC 
analyzed fields from 1200Z 10 May 1973. 

As indicated in section III the angular phase speed for wave number 
six was “9® longitude/day (toward the west). A 72-hour forecast of 
surface pressure was made by both the second and mixed order differenc- 
ing. The changes in phase angle vs latitude were plotted for both 

f 

cases and are shown in Figures 3 and 4. Table I summarizes the average 
angular movement per 24-hour the entire 72-hour forecast for both dif- 
ferencing schemes. It should be noted that the mixed differencing 
scheme propagated the wave slower to the west than the second order 
scheme. The barotropic Rossby equation (Holton, 1972) is given as 

c - S - ■ / 

(k^+ m^) 



where B is , f is the Coriolis parameter, u is the mean zonal flow, 
k is the zonal wave length, m is the wave number and C is the zonal 
phase speed. Since the quantities B> k and m are constant the phase 
speed, C, depends only upon the mean zonal flow. Fourth-order differ- 
encing should give a better approximation for the effect of the mean 
zonal wind, u, and therefore the westward movement of the wave would 
be slower than that of the second order differencing case. It would 
appear that the mixed differencing scheme does give a better solution. 
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Figure 3. Phase Angle vs Latitude for the Second Order Scheme. 
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Figure A. Phase Angle vs Latitude for the Mixed Scheme. 
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LATITUDE 


SECOND ORDER 
DEG LAT/DAY 


MIXED SCHEME 
DEG LAT/DAY 


0 

O 


-12.6® 


-4.8 


60® 


-10.3® 


-7.9 


50® 


O 

CM 

• 

o\ 

1 


-7.4 


O 

O 


-8.5® 


-7.5 


30® 


0 

• 

1 


-7.2 


20® 


-7.3® 


-7.6 


10® 


-7.0® 


-6.5 


0® 


-4.6® 


-3.6 



TABLE I. Average Movement of Analytic Wave per 24-Hours for a 

72-Hour Forecast for Both Second Order Differencing and 
the Mixed Scheme. 
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A 48-hour forecast was also made by both the second order and 
mixed differencing schemes using analyzed FNWC data fields. Chart A 
is the initial analyzed surface pressure field. Charts B through E 
are the 12- , 24-, 36- and 48-hour forecasts produced by the second 
order scheme and Charts F through I are the corresponding forecasts 
produced by the mixed differencing scheme. Table II lists five high 
pressure centers and Table III lists seven low pressure centers with 
their initial position and both the 48-hour second order and 48-hour 
mixed scheme forecast positions. As can be seen each high pressure 
center was forecast at least five degrees farther to the east by the 
mixed differencing scheme. Low pressure centers 1 and 3 also were 
forecast five degrees farther to the east by the mixed scheme while 
low centers 2 and 4 were moved very little by both differencing schemes. 
Low number 5 was the only low that was moved slower than the second 
order scheme. Low number 6 was developed into two centers by both 
forecast schemes. Only the forecast position of the major center was 
listed in Table III. The mixed scheme moved this center nine degrees 
farther than the second order scheme. Center number 7 was also 
initially a single center but was developed into three small centers 
with the second order forecast scheme. The mixed scheme maintained 
the more realistic single center system near its initial position. 
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Chart C. 



24-Hour Forecast, 



Second Order Scheme, 
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Chart D. 36-Hour Forecast, Second Order Scheme. 
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Chart E. 48-Hour Forecast, Second Order Scheme. 
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Chart H. 



36-Hour Forecast, 



Mixed Scheme. 
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Chart I. 



A8-Hour Forecast, Mixed Scheme. 
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HIGH CENTER 
NUMBER 


INITIAL 

POSITION 


48-HOUR SECOND 
ORDER POSITION 


48-HOUR MIXED 
SCHEME POSITION 


• 1 


37N 177W 


37N 17 7W 


35N 


172W 


2 


35N 136W 


47N 128W 


43N 


12 IW 


3 


None 


31N 53W 


33N 


48W 


4 


45N 


25W 


52N 17W 


52N 


12W 


5 


64N 


50E 


60N 60E 


52N 


66W 



TABLE II, Initial Positions of High Centers and Forecast Positions 
of the Second Order Scheme and the Mixed Scheme. 



LOW CENTER 
NUMBER 


INITIAL 

POSITION 


48-HOUR SECOND 
ORDER POSITION 


48-HOUR MIXED 
SCHEME POSITION 


1 


47N 155W 


48N 152W 


49N 


14 7W 


2 


47N 


89W 


SON 4SW 


48N 


40W 


3 


47N 


SOW 


SON 87W 


SIN 


88W 


4 


67N 


02E 


68N 03E 


68N 


04E 


5 


61N 


89E 


S9N 99E 


61N 


97E 


6 


36N 152E 


SON ISIE 


48N 


160E 


7 


30N 


80E 


* 


33N 


81E 



TABLE III. Initial Positions of Low Centers and Forecast Positions 
of the Second Order Scheme and the Mixed Scheme. 



^Several small lows were forecast by the second order scheme. 
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VI. CONCLUSIONS 



Both the analytic and real data experiments indicate that the 
mixed fourth and second order differencing scheme gives a more 
accurate solution for the meteorological phase speeds. More analytic 
experiments with a more realistic phase speed will have to be per- 
formed to produce the final conclusive proof of the mixed scheme. 

A series of real data experiments should also be considered. In 
addition the fourth order differencing must eventually be extended 
to the poles. 
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APPENDIX A 



DIFFERENCE EQUATIONS 



The longitudinal momentum equation terms in difference form are 



n+1 n-1 

dUTT _ l.i IJ 

9t “ 2At ’ 



( 1 ) 



1 r9(uU7T) 9 (uVTT)cOS0, _ 
a COS 0 ^ 9X 90 ^ ~ 



. (u..+ u.,- .)U1T.., /„ . - (u..+ u. 1 .)u7T. , . 

1 ■> 11 1+1, r ij 1-1, 1-1/2, 1 

a cos 0. ^ 2AX 






where 



U1T.,,,_ . = -^(iT. .U. .+ 77. , , -u.,, .+ 77. ..,U. .,,+ 77. . ,U. . ,) , 

X+1/2,J 4 XJ xj X+1,J X+1,J 1,J+1 x,j+l X,J-1 X,J-1 

( 2 . 1 ) 

U7T. , . = - 7 ( 77 . .U. .+ 77. , .U. , .+ 77. , .,,U. , , 

1-1/2, j 4 XJ XJ X-1,J X-1,J X-1,J+1 x-l,j+l + 



i. , . ,u. , . ,) 

x-l,j-l x-l,j-l 



( 2 . 2 ) 



U7T 



i. j+1 



= - 7 ( 77 . .U. .+ 77. ..„U. ...,+ 77. , .,iU. , + 

4 IJ IJ i»J+2 i,J+2 x-l,j+l x-l,j+l 



77. . ,U. . ,) 

1,J-1 l.J-l 



(2.3) 
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(2.4) 



UTT. . , = •tC'JT. .U..+ IT. . -U, . -+ TT. , . ,U . 1 + 

l,j-l 4' ij IJ i,j-2 i,j-2 



IT. . tU . t) 

i,3-l 1,3-1 



The longitudinal grid increment (AX) and the latitudinal grid increment 
(A0) are both ten degrees. 



- W-I M ('T. • 1 . 1 + u. .1 ) w. . , , (u. . + u. . ) 

9wu _ - , ink i3,k+l 13k _ X3,k-1 13k xi,k-l 

9a ■ ij ^ 2Aa 2Aa 



(3) 



where the vertical velocity, -O , is defined as. 



w, 



k+1 



Aa ,3 it,N 

Ij 



, UTT.,, - UTT. -UTT... VTT. . , COS0. , ,-vi. . ,,COS0. , 

1 , 1+1, lk x+l.ik xik xi+lk 1+1 xi-lk i~l 

a cos 0. ^ AX A0 

^ (4) 



and the local change as 



(- 



UTT -UTT. . 

r i+llk X3k 

'i-j / j AXa cos 0. 



_ 

at^ij 



k=l 



VTT. . COS0 . , ,-VTT. . ,,COS0. , 
X3+lk 3+1 X3-Ik 3-1 . 

A0 a cos 0 . ^ 

J 



(5) 



The UTT terms of Equations (4) and (5) are. 



UTT. . = U. . , . (tT. .+ TT. , , TT, , ,,+ 
X+I3 1+I3' X3 X+I3 X3+I 



TT. . ) -i 

13-1 4 



(5.1) 



The quantities titt^^ , and are also determined in the 



same manner. 
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■A. 






( 6 ) 



wf = fj 



where 



"ij ■ 



TT. , .+ TT. T . ,-+ TT. , . ,) 

1-lj 1-1, J+1 1-1, J-1 



( 6 . 1 ) 



TTuv tan9 
a 



u . . (ttv . . ) tan9 . 

J-3 1 

a 



(7) 



1 rRT 3tT ^ TT 
a COS e ^ 3X 3X^ 



a cos 



0 . 

J 



fRT (— ) + TT — iziiJ—] 

^ ij ’'3X-' ij AX ^ 



( 8 ) 



where 



T.. E -7(1. T.+ T..+ T. , .^1+ T. , . t) 

ij ^ i-lj iJ 1-1, J+1 1-1, J-1 



i . . = -r(Tr . ,.+ 7r..+7r. , , . ,) 

ij 4' i-lj ij 1-1, J+1 1-1, J-1' 



( 8 . 1 ) 

( 8 . 2 ) 



From Equation (9) of section II the change of terrain height is found 
as follows 



(il )* = IL_ r. 
''3X^^ RT '■ 



= — [■ 

RT '■ 











AX 


P “ 


AX 




- 




(, . 




E_ 


a 


0 





AX 



a] 




(9) 

( 10 ) 



i 
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where $ is to the right of the computational point and $ is to the 

R. li 

left of the computational point. Sigma and p indicate the surface in 
question. The thickness between the pressure surface and the sigma 
surface to the right and left of the computational point is defined as. 





and 



( 11 ) 





( 12 ) 



therefore one can write. 







( 13 ) 



( 14 ) 



The mean temperature is found as follows. 



T = 



T 

comp 



3T,A In 7T, 
30^ 2 ^ 



( 15 ) 



and in difference form 





(T 



“ T ) 

R I. R 
0+1 0 

A0 




(T 



^ 0+1 

Aa 



(LN IT . . - LN TT . . ) 

ii_ 

2 



(LN IT . . - LN 7T . . ) 

. il_ 

2 



( 16 ) 

( 17 ) 
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(17.1) 



T = — fT + T + T + T ) 

ij ±3 i-lj i-lJ+l 



LN 7T , . = 7T, . + Jin tt. - . + Jin 7T , - . , - + Jin tt. - . - ) (17.2) 

ij ^ ij i“l3 



Figure 5 illustrates the differencing in both the horizontal and verti- 
cal* Substituting Equations (16) and (17) into Equations (13) and (14) 
one obtains. 



^ (Hn TT,.-Jln Tr..) 

"r S ' —^^2 ) (18) 

op '* J J 



= -R[T. .4 
L ij 

op 



/■t _T ■) _ 

'■ L .1 (£.n IT.. -Jin it. , .) 

0+1 o n i-li 



Aa 



] Jlmr..-Jlmr. .. (19) 

ij 1 -I 3 



where 






= Jin 



0+1 



(°+l . - 

= Jin TT. . - Jin TT. . 



( 20 ) 



and 



Jln_ = Jin TT . . - Jin TT . ^ . 

L P^+^ ij 1 -I 3 



( 21 ) 
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Figure 5. Interpolation for Vtt. 
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The latitudinal-momentum equation terms in difference form are 



N+1 N-1 

~ TTv, . - irv, . 

3vtt _ ij ij 



3t 



2At 



(1) 



1 ^ 3 (ttvw) ^ 3 (ttwcosG) j 



a cos 0 3X 



30 



, (v. .+V. ,, .)uTT. , , - (v..+v. T .)u7T. . 

1 r 13 • 1+1-1 i+l/2i n 1-l.r 1-1/2.1 



a cos 0 . 

3 



2AX 



(v..+v..,„)vir. COS0 - (v..-v.. ^)v7T.. cos0. 

. 13 13+2 13+1 3+1 13 13-2 13 ~1 Jzli (o) 

2A0 ^ ^ ^ 



The uTT terms are determined in the same manner as the terms for 
Equation (2) were found. 



7r8(wv) _ - r g iia+1 ijg g~l 13 a ijg-l ^ 

3o " ”^ij^ 2Aa ~ 2Aa ^ 



(3) 



The vertical velocity, w, and the local change, (“^) . . 9 

dt ij 

as in Equations (4) and (5) of the longitudinal equation 



are obtained 
of momentum. 



irvf = f . TT . . V . . 
1 ij ij 



(4) 



TT. . = Equation (6.1) or the u-equation of momentum 4* 



7TUU tan 6 
a 



u. . (tt. .u. .) tan 0 . 

1 



(5) 



V 



46 



i fTi ^ . RT _3rr, 

“ a ^ 36 30^ 



a cos 0j ij 39^ ij A0 



- V> <V - *L > 

,3tTv* _ Ti_ r p g _a £_ 

RT ^ A0 



] 



( 6 ) 



(7) 



Equation (7) is solved by the same method as was shown for the 
lontitudinal momentum equation. 

Thennodynamic Equation 



9(ttT) 3 ^’^ija ~ '^’^ija 
3t ~ 2At 



( 1 ) 



1 ,3(ttuT) 3(ttvTcos9 , 
I J 



a cos 6 3A 



30 



, T. , T . (u. , - .TT_ , .) - T..(u..tt..) 

1 . l+lj 1+1.1 1+1 J IJ 3-3 IJ 

a cos 6 . AX 



)cos0 +T.. tV.. it.. cos0. 
1,1+1^ ij+1 x.1+1 j+1 ij-l ij-1 ij-1 j-1 

A0 



1 (2) 



where T is the temperature at the velocity point based on the average 
of the four surrounding mass points and it is also calculated in the 
same manner. 
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TT 3(wT) _ _ r ”l1o^\io+l “l1o-l^ 

do ~ ^ij ^ 2Ao ~ 2Ao 



] (3) 



The vertical velocity, w, and the local change, (“^)y » obtained 

as in the u-equation of motion. 



RT r ^ 1 r ^ 

r -WTT+aC;^ + + V COS0J) 

c 0 I dt a cos 9 dA 39 



]■ 



(w + W T ) 
r' O 0-r „ . „ 

c„0 ^2 ij a 



RT 



(— ) 



+ 



0 

0 



a cos 9. 

3 



(u, , - . 77. , - . - U. . TT. .) 

r i+lj i+lj iJ LL_ 

^ 2AA 



+ 



(v. . . ^TT. . .^COsG. -V.. ,TT.. ^COS0. t) 

11+1 11+1 1+1 1 , 1-1 n -1 , 1-1 . 

2A0 ^ 



(4) 



where the vertical velocity and the local change are again computed 
as before. 



The Moisture Equation 



N+1 N-1 

„ qTT . . - qTT . . 

3qTT _ ij 11 

3t “ 2At 



( 1 ) 



1 , 3uTTq 3 (vTfqcosG) .. _ 

a cos 0 3A 30 ~ 



a cos 0 . 

J 



q . . (u. .ir_ , .) - q..(u..Tr..) 

^i+lj i+lj i+lj ^11 ij ij 

AX 



qj --LI (v ... .IT. )cos0 - - q.. ,(v.. ,ir.. ,cos0. ,) 
ll+l 11+1 11+1 1+1 ^11-1 ii-l 3-1-1 1-1 

Ao 



(2) 
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V 



■ ■ 






# 




t 



where q is the specific humidity at the velocity point based on the 
average of the four surrounding mass points and tt is also calculated 
in the same manner. 



IT 3(wq) _ p ^iia^^ijg+l 

do ij '■ 2Aa 



w.. ,(u.. +u.. ,) 

ng-1^ no 

2Aa ^ 



(3) 



The Hydrostatic Equation 



d(I) = - RT diln O 



6 = $ - R£n [ — n — n ] 

k+1 k ^ ^^^k+1 ~ 2£n7T 



At the lowest level $ is given as 



. = 7T - RT. £n (0.9) 
ij ijk 
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APPENDIX B 



In accordance with the work of Williams (1972) the fourth order 
terms are a sum of a weighted average of the second order and fourth 
order differencing. The weighting factor the second order part is 
four-thirds and a minus one-third for the fourth order portion. The 
advection term of the longitudinal momentum equation is: 

- / (u. .+ U. , - .)uTT. - (u..+ u. - .)U7T. - . 

1 r n 1+1.1 1 - 1/23 1-lj l-l/2j 

a cos 0. 4 ^ 2AX 

J 



(u..+ u. . , ^)cos 0 . , - VTT..,. -(u..- U.. «)cos 0 . - V 7 T . . 

+ 3-3+2 1+1 ij+1 ij ij-2' iJzl) 

2AX ^ 



, (u. .+ U.j „ .)uTT.j, .-(u. .+ U. „.)UTT. 

Ir ^ 11 l+2j' l+ll 11 X-2.1 1-1,1 

~ 3^ 4AX 



(u..+ u. . , . )cos9. VTT. -(u..+ u.. ,)cos9. „ VTT.. _ 

I 13 13+^ 3+2 13+2 IJ ±2^ m 

4AX ^ 



The VTT terms are determined as shown in Appendix A. The advection term 
for the latitudinal momentum equation is similar to term (1). The 
advection term of the thermodynamic equation is: 

i r l+lj l+lj 13 IJ 



[t { 



AX 



a cos 0. 3 

J 

T. . , ,U7T. . COS0 . , ^ -T.. - UTT.. COS0. - 
13+1 lJ+1 J+1 lJ-1 lJ-1 J-1 t 

AX ^ 

“IT* 

1 T. , , . UTT . - T. , . uir. . . 

_ Ir l+lj i+lj izll + 

3^ 2AX 



— ^ ^ —it 

T UTT. cos9 - T.. „ uir.. cos9 . 

13+2 13+2 j+2 13-2 13-2 j-2 

2A9 



}] 



( 2 ) 



V 



50 







i 



The program as written gives the average temperatures, T, and utt values 



at velocity points for the second order differencing. For fourth order 
differencing the grid distance is doubled and the values T and utt are 



determined at the adjacent mass points. All the averaged quantities 



are determined similarly to those in Appendix A. 



The flux term of the continuity equation is: 



-I / • I 1 • 1 “ UTT . . , 

1 A f 1+1.1 k Ilk 

a cos 9.3 AX 

3 



vir. , , T, COS0 . -vTT.. ^,cos9. 
xi+lk 1+1 1 , 1 -lk i-l -, 

AX ^ 



~T* “=■* 

T UTT. - UTT. 

i f 1+1.1 k 1-llk 
" 3 ^ 2AX 



— 

VTT. . , COS0 . , „ - VTT.. -COS0. „ 
ij+2k .1+2 1-1-2 .1-2 

2A0 



}] 



The averaged values are also determined as indicated in Appendix A. 
The unstarred quantities are at a velocity point and the starred 
quantities are at a mass point. For this experiment the fourth order 
differencing was blended with the second order differencing at 75® 
latitude. This was done to avoid having to go over the poles to pick 
up the j+1, j+2 and j+3 quantities at the higher latitudes. 
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